
rm(list=ls(all=TRUE))
library(ggplot2)
library(gridExtra)
library(car)
library(plm)
library(gam)
library(plyr)
library(doBy)
library(foreign)
library(texreg)
library(xtable)
library(MASS)
library(GGally)
library(zoo)
library(lmtest)
library(plm)
library(lattice)
library(ebal)
library(DataCombine)
library(rgenoud)
library(sandwich)
library(Hmisc)
library(dplyr)
library(tidyr)
library(lubridate)
library(zoo)
require(texreg)
require(lmtest)
require(multiwayvcov)
library(VIM)
library(rpart)

setwd("/users/irismalone/Dropbox/AGDPapers/IS-ISSS Paper/Replication Files - II/Data")

df2= read.csv("replication_extsup_data.csv")
dim(df2)

########################
# Summary Statistics
########################


df = df2 %>% distinct(ccode1, ccode2, torgid, .keep_all = TRUE)

length(unique(df$torgid))
length(unique(df$torgid[df$anystatesup==1]))
length(unique(df$torgid[df$civwar==1]))
length(unique(df$ccode))
length(unique(df$ccode2[df$anystatesup==1]))

require(stargazer)
df$sharedcoideo = df$leftideational+df$sharedrelig+df$sharedethnicpower
df$sharedcoideo[df$sharedcoideo >0] = 1

df = subset(df, df$torgid!=28 & df$polopp==1)#remove Al Qaeda outlier and subset to polopp

dfsummstats = with(df, data.frame(anystatesup, sharedcoideo, sharedethnicpower, leftideational, sharedrelig, shiaideational, sunniideational))
stargazer(dfsummstats, median = T, omit.summary.stat = c("p25", "p75"))


########################
# ROBUSTNESS CHECKS
########################


######
# Table 3
# Material Support
######

df3=df2
df3$anystatesup = ifelse(df2$anystatesup==1 & (df2$statetypesupmaterial== 1 | 
                                                 df2$statetypesuptrain== 1 | 
                                                 df2$statetypesupfinancial== 1 | 
                                                 df2$statetypesupother== 1), 1, 0)

reg.m0 = glm((anystatesup) ~ l1noviolsupl + sharedcoideo, 
             data=df3, family="binomial")
vcov_firm <- cluster.vcov(reg.m0, df3$dyadid)
a = coeftest(reg.m0, vcov_firm)
a

reg.m1 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo, 
             data=df3, family="binomial")
vcov_firm <- cluster.vcov(reg.m1, df3$dyadid)
b = coeftest(reg.m1, vcov_firm)
b

reg.m2 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo  + rivalry + minwbgdppc + mindemest, 
             data=df3, family="binomial")
vcov_firm <- cluster.vcov(reg.m2, df3$dyadid)
c = coeftest(reg.m2, vcov_firm)
c


reg.m3 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo + rivalry + minwbgdppc + mindemest + factor(ccode1), 
             data=df3, family="binomial")
vcov_firm <- cluster.vcov(reg.m3, df3$dyadid)
d = coeftest(reg.m3, vcov_firm)
d

texreg(list(a, b, c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Material Support Only")
texreg(list(reg.m0, reg.m1, reg.m2, reg.m3))

htmlreg(list(a, b, c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Material Support Only", file = "table3_appendix_materialsupp.doc")

htmlreg(list(reg.m0, reg.m1, reg.m2, reg.m3), file = "table3_appendix_materialsupp_unclustered.doc")


######
# Table A4
######

####################
# Substitutability - Same Movement
####################

reg.m0 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo2 + sharedethnicpower  + rivalry + minwbgdppc + mindemest, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m0, df2$dyadid)
a = coeftest(reg.m0, vcov_firm)
a

reg.m1 = glm((anystatesup) ~ l1noviolsupl*sharedethnicpower + sharedcoideo2  + rivalry + minwbgdppc + mindemest, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m1, df2$dyadid)
b = coeftest(reg.m1, vcov_firm)
b

reg.m2 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo2 + sharedethnicpower  + rivalry + minwbgdppc + mindemest + factor(ccode1), data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m2, df2$dyadid)
c = coeftest(reg.m2, vcov_firm)
c


reg.m3 = glm((anystatesup) ~ l1noviolsupl*sharedethnicpower + sharedcoideo2 + rivalry + minwbgdppc + mindemest + factor(ccode1), data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m3, df2$dyadid)
d = coeftest(reg.m3, vcov_firm)
d


texreg(list(a, b, c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Substitution")
texreg(list(reg.m0, reg.m1, reg.m2, reg.m3))

htmlreg(list(a, b, c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Different Types of Ties", file="table4_appendix.doc")
htmlreg(list(reg.m0, reg.m1, reg.m2, reg.m3), file="table4_appendix_unclust.doc")

####################
# Table A5
# Substitutability - Rebel Groups
####################

reg.m0 = glm((anystatesup) ~ l1noviolsupl + civwar + sharedcoideo, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m0, df2$dyadid)
a = coeftest(reg.m0, vcov_firm)
a

reg.m1 = glm((anystatesup) ~ l1noviolsupl*civwar + sharedcoideo, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m1, df2$dyadid)
b = coeftest(reg.m1, vcov_firm)
b

reg.m2 = glm((anystatesup) ~ l1noviolsupl*civwar + sharedcoideo  + rivalry + minwbgdppc + mindemest, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m2, df2$dyadid)
c = coeftest(reg.m2, vcov_firm)
c

reg.m3 = glm((anystatesup) ~ l1noviolsupl*civwar + sharedcoideo + rivalry + minwbgdppc + mindemest + factor(ccode1), data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m3, df2$dyadid)
d = coeftest(reg.m3, vcov_firm)
d


texreg(list(a, b, c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Alt Substitution Measures")
texreg(list(reg.m0, reg.m1, reg.m2, reg.m3), stars=c(0.1, 0.05, 0.01))

htmlreg(list(a, b, c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Alt Substitution Measures", file="table5_appendix.doc")
htmlreg(list(reg.m0, reg.m1, reg.m2, reg.m3), file="table5_appendix_unclust.doc")



####################
# Alternative Country Controls
####################

reg.m2 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo  + l1nostatesupl + rivalry + cowmaj2 + 
               mindemest + minwbgdppc + lndist , data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m2, df2$dyadid)
c = coeftest(reg.m2, vcov_firm)
c


reg.m3 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo  + l1nostatesupl + rivalry + cowmaj2 + 
               mindemest + minwbgdppc + lndist + factor(ccode1), data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m3, df2$dyadid)
d = coeftest(reg.m3, vcov_firm)
d

texreg(list(c,d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Additional Country Controls")
texreg(list(reg.m2, reg.m3), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01))

htmlreg(list(c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Alt Country Controls", file="table6_appendix.doc")
htmlreg(list(reg.m2, reg.m3), file="table6_appendix_unclust.doc")

####################
# Alternative Group Controls
####################

reg.m2 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo  +  rivalry  + 
               mindemest + minwbgdppc + aim2 + polwing + transnationalattackany + civwar, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m2, df2$dyadid)
c = coeftest(reg.m2, vcov_firm)
c


reg.m3 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo  +  rivalry  + 
               mindemest + minwbgdppc + aim2 + polwing + transnationalattackany + civwar + factor(ccode1), data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m3, df2$dyadid)
d = coeftest(reg.m3, vcov_firm)
d

texreg(list(c,d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Additional Group Controls")
texreg(list(reg.m2, reg.m3), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01))


htmlreg(list(c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Alt Group Controls", file="table7_appendix.doc")
htmlreg(list(reg.m2, reg.m3), file="table7_appendix_unclust.doc")


########################
# Modeling Specifications
########################

#Cross-Sectional FE

reg.m2 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo  + rivalry + minwbgdppc + mindemest + factor(ccode2), data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m2, df2$dyadid)
c = coeftest(reg.m2, vcov_firm)
c


reg.m3 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo + rivalry + minwbgdppc + mindemest + factor(region.x), data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m3, df2$dyadid)
d = coeftest(reg.m3, vcov_firm)
d

texreg(list(c, d), omit.coef="ccode2", stars=c(0.1, 0.05, 0.01), caption="Alternate Panel Fixed Effects")
texreg(list(reg.m2, reg.m3))

htmlreg(list(c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Alt Panel Fixed Effects", file="table8_appendix.doc")
htmlreg(list(reg.m2, reg.m3), file="table8_appendix_unclust.doc")


##############
# Table 9
# Control for IS
##############
table(df2$coldwar)
reg.m0 = glm((anystatesup) ~ l1noviolsupl + sharedcoideo + coldwar + post911, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m0, df2$dyadid)
a = coeftest(reg.m0, vcov_firm)
a

reg.m1 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo + coldwar + post911, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m1, df2$dyadid)
b = coeftest(reg.m1, vcov_firm)
b

reg.m2 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo  + rivalry + minwbgdppc + mindemest + coldwar + post911, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m2, df2$dyadid)
c = coeftest(reg.m2, vcov_firm)
c

reg.m3 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo + rivalry + minwbgdppc + mindemest + coldwar + post911 + factor(ccode1), data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m3, df2$dyadid)
d = coeftest(reg.m3, vcov_firm)
d


texreg(list(a, b, c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Intl System Controls")
texreg(list(reg.m0, reg.m1, reg.m2, reg.m3), stars=c(0.1, 0.05, 0.01))

htmlreg(list(a, b, c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="International System Periods", file="table9_appendix.doc")
htmlreg(list(reg.m0, reg.m1, reg.m2, reg.m3), file="table9_appendix_unclust.doc")


#######
# Table 10
# Linear Time Trend
#######
reg.m0 = glm((anystatesup) ~ l1noviolsupl + sharedcoideo + year, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m0, df2$dyadid)
a = coeftest(reg.m0, vcov_firm)
a

reg.m1 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo + year, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m1, df2$dyadid)
b = coeftest(reg.m1, vcov_firm)
b

reg.m2 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo  + rivalry + minwbgdppc + mindemest + year, data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m2, df2$dyadid)
c = coeftest(reg.m2, vcov_firm)
c

reg.m3 = glm((anystatesup) ~ l1noviolsupl*sharedcoideo + rivalry + minwbgdppc + mindemest + year*factor(ccode1), data=df2, family="binomial")
vcov_firm <- cluster.vcov(reg.m3, df2$dyadid)
d = coeftest(reg.m3, vcov_firm)
d


texreg(list(a, b, c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Linear Time Trend")
texreg(list(reg.m0, reg.m1, reg.m2, reg.m3), stars=c(0.1, 0.05, 0.01))

htmlreg(list(a, b, c, d), omit.coef="ccode1", stars=c(0.1, 0.05, 0.01), caption="Linear Time Trend", file="table10_appendix.doc")
htmlreg(list(reg.m0, reg.m1, reg.m2, reg.m3), file="table10_appendix_unclust.doc")

#######
# Table 11
# Rare Event Logit
#######

require(Zelig)
z.out1 =  zelig(anystatesup ~ l1noviolsupl + sharedcoideo,
                data = df2, model = "relogit", tau = 1042/303772)

z.out2 =  zelig(anystatesup ~ l1noviolsupl*sharedcoideo,
                data = df2, model = "relogit", tau = 1042/303772)

z.out3 =  zelig(anystatesup ~ l1noviolsupl*sharedcoideo  + rivalry + minwbgdppc + mindemest,
                data = df2, model = "relogit", tau = 1042/303772)

z.out4 = zelig(anystatesup ~ l1noviolsupl*sharedcoideo  + rivalry + minwbgdppc + mindemest + factor(ccode1),
                data = df2, model = "relogit", tau = 1042/303772)

texreg(list(z.out1, z.out2, z.out3, z.out4), omit.coef="ccode1")
htmlreg(list(z.out1, z.out2, z.out3, z.out4), file="table11_appendix.doc")


# Sampling
require(caret)
set.seed(2022)
cvseeds = vector(mode = "list", length = 56) 
for(i in 1:55) cvseeds[[i]] <- sample.int(1000, 56) 
cvseeds[[56]] <- sample.int(1000, 1) 
train_control_smote = trainControl(method="repeatedcv", #method
                             number=10, #number of folds
                             repeats = 5,
                             savePredictions=T, 
                             classProbs=T,  # should class probabilities be returned
                             sampling="smote",  #smote sample
                             summaryFunction=twoClassSummary, 
                             seeds=cvseeds)

df2$statesup=as.factor(df2$anystatesup)
levels(df2$statesup)=c("Yes","No")

glm.model.smote = train(statesup ~ l1noviolsupl*sharedcoideo  + rivalry + minwbgdppc + mindemest,
                        data=df2,
                        method = "glm", 
                        trControl=train_control_smote,
                        na.action=na.omit)

train_control_up = trainControl(method="repeatedcv", #method
                                   number=10, #number of folds
                                   repeats = 5,
                                   savePredictions=T, 
                                   classProbs=T,  # should class probabilities be returned
                                   sampling="up", 
                                   summaryFunction=twoClassSummary, 
                                   seeds=cvseeds)

glm.model.up = train(statesup ~ l1noviolsupl*sharedcoideo  + rivalry + minwbgdppc + mindemest,
                        data=df2,
                        method = "glm", 
                        trControl=train_control_up,
                        na.action=na.omit)


train_control_down = trainControl(method="repeatedcv", #method
                                number=10, #number of folds
                                repeats = 5,
                                savePredictions=T, 
                                classProbs=T,  # should class probabilities be returned
                                sampling="down", 
                                summaryFunction=twoClassSummary, 
                                seeds=cvseeds)

glm.model.down = train(statesup ~ l1noviolsupl*sharedcoideo  + rivalry + minwbgdppc + mindemest,
                        data=df2,
                        method = "glm", 
                        trControl=train_control_down,
                        na.action=na.omit)

texreg(list(glm.model.down$finalModel, glm.model.up$finalModel, glm.model.smote$finalModel), stars=c(0.1, 0.05, 0.01))
htmlreg(list(glm.model.down$finalModel, glm.model.up$finalModel, glm.model.smote$finalModel), file="table12_appendix.doc")

